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A  nonlinear  finite  element  procedure  has  been  developed  for  the  dynamic 
vibrational  analysis  of  planar  mechanisms.  The  analysis  takes  into  account  the 


effects  of  geometric  and  material  nonlinearities,  vibrational  effects  and 
coupling  of  deformations.  Numerical  results  have  been  reported  for  certain 
mechanism  examples.  The  effects  of  the  nonlinearities  have  been  found  to  be 
significant  on  the  dynamic  behavior.  Due  to  the  complex  nature  of  this 
nonlinear  analysis  procedure,  an  efficient  optimal  design  approach  using  an 
optimality  criterion  technique  has  also  been  developed.  The  new  optimization 
technique,  called  the  Gauss  Nonlinearly  Constrained  Technique,  has  been 
developed  in  such  a  way  that  it  is  applicable  to  design  problems  with  nonlinear 
objective  functions  and  constraints.  The  applicability  of  this  method  has  been 
demonstrated  with  example  problems  consisting  of  objective  functions  of  various 
complexities.  Complete  details  of  the  nonlinear  finite  element  procedure  as 
well  as  the  optimization  technique  are  available  in  two  papers  (to  be  published) 
which  are  included  here  in  the  Appendix. 


The  objectives  of  this  research  are  to  develop  a  nonlinear  finite  element 
dynamic  analysis  procedure  for  both  planar  and  spatial  mechanisms  of  the  type 
found  in  space  structures.  Geometric,  and  material  nonlinearities  as  well  as 
combinations  of  these  are  to  be  included  in  the  analysis.  Furthermore,  an 
efficient  optimal  design  method  is  to  be  developed  to  handle  objective  functions 
composed  of  a  combination  of  rigid  body  and  deformation  displacements  involving 


geometric  design  variables  es  veil  as  cross-sectional  sizes.  Thus,  during  the 
reporting  period  it  vas  proposed  to 


a)  develop  a  nonlinear  finite  element  procedure  for  dynamic  planar 
mechanisms; 

b)  develop  an  efficient  optimization  method  involving  a  small  number  of 
analyses  for  mechanism  design  problems. 


Substantial  progress  has  been  made  in  both  of  these  two  research  areas 
during  the  reporting  period.  Complete  details  of  the  work  are  included  in  the 
Appendix;  only  a  summary  is  presented  here. 

In  the  nonlinear  finite  element  analysis  area,  planar  mechanism  problems 
have  been  formulated  to  include  the  effects  of  material  and  geometric 
nonlinearities  due  to  large  deformations,  and  vibrational  effects  for  members 
with  nonuniform  cross  sections.  Effects  of  transverse  shear  and  rotatory 
inertia  have  also  been  included.  Numerical  results  for  deformationo  and 
stresses  have  been  reported  for  certain  planar  mechanisms.  The  results  clearly 
indicate  the  effects  of  geometric  and  material  nonlinearities  on  the  dynamic 
behavior  of  planar  mechanisms  and  the  need  to  include  these  effects  in  the 
analysis. 

In  the  optimal  design  area,  a  new  optimization  technique,  called  the  Gauss 
Nonlinearly  Constrained  Method  (GNLC),  which  is  applicable  to  design  problems 
with  nonlinear  objective  functions  and  constraints  has  been  developed.  The  GNLC 
technique  is  an  extension  of  a  previously  developed  method  called  the  Gauss 
Constrained  Technique  which  is  capable  of  handling  only  linear  constraints. 

Both  of  these  techniques,  based  on  Gauss'  unconstrained  method,  were  developed 


so  that  the  Kuhn-Tucker  conditions  were  automatically  satisfied  when  the 
procedure  terminated.  As  an  example,  a  four-bar  mechanism  was  designed  such 
that  the  coupler  point  (XX,YY)  of  the  mechanism  closely  generated  a  curve 
defined  by  discrete  points.  The  objective  function  to  be  minimized  in  this  case 
was  the  sum  of  the  distances  (squared)  between  the  desired  curve  and  the  actual 
curve  generated  by  the  mechanism.  This  objective  function  is  a  highly  nonlinear 
function  of  the  design  variables.  For  the  particular  example  considered,  the 
GNLC  technique  was  found  to  require  significantly  fewer  objective  function  and 
constraint  evaluations  than  currently  available  methods. 
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THE  DEVELOPMENT  AND  APPLICATION  OF 
GAUSS'  NONLINE ARLY  CONSTRAINED  OPTIMIZATION  METHOD 

D.R.  Boston*,  K.D.  Willmert* 
and  M.  Sathyamoorthy* 


Abstract 

Presented  in  this  paper  is  a  new  optimization  technique,  called 
the  Gauss  Nonlinearly  Constrained  Technique  which  is  applicable  to  de¬ 
sign  problems  with  nonlinear  objective  functions  and  constraints.  The 
technique  is  an  extension  of  a  previously  developed  method  for  linear 
constraints,  referred  to  as  the  Gauss  Constrained  Technique.  Both  of 
these  techniques,  based  on  Gauss'  unconstrained  method,  have  been  de¬ 
veloped  so  that  the  Kuhn-Tucker  conditions  are  automatically  satisfied 
when  the  procedure  terminates. 

Introduction 

The  optimal  design  of  many  structures  and  mechanical  mechanisms  in¬ 
volves  one  or  more  complex  and  time  consuming  analyses  at  each  iteration 
of  the  optimization.  This  is  particularly  critical  if  a  large  deformation 
nonlinear  analysis  is  required.  In  these  cases  especially,  it  is  im¬ 
portant  that  the  optimization  method  require  very  few  analyses,  even  at 

the  expense  of  significantly  increasing  the  amount  of  calculations  by  the 
optimization  technique  itself.  For  example  in  a  recent  work,  DeRubes 
and  Willmert  [1]  applied  the  relatively  efficient  Generalized  Reduced 
Gradient  (GRG)  Technique  of  Lasdon  et.al  [2]  to  mechanism  design  for  path 
generation  and  rigid  body  guidance.  The  mechanism  links  were  considered 

♦Mechanical  and  Industrial  Engineering  Department,  Clarkson  University, 
Potsdam,  N.Y.  13676. 


flexible,  and  thus  a  quasi-static  (linear)  finite  element  analysis  was 
used  to  obtain  deformations  and  stresses.  The  GRG  required  as  many  as 
14,000  mechanism  analyses  to  obtain  the  optimal  design.  Computation 
time  approached  twenty  hours  on  an  IBM  4341  mainframe  computer.  If  a 
nonlinear  analysis  had  been  used,  the  corresponding  times  would  have 
been  considerably  higher. 

To  reduce  the  number  of  analyses,  Paradis  and  Willmert  [3]  develop¬ 
ed  a  new  direct  method  for  efficient  design  of  mechanisms.  Gauss' 
method,  which  Wilde  [4]  concluded  to  be  very  efficient  for  unconstrain¬ 
ed  mechanism  design,  was  modified  to  handle  linear  constraints.  The  re¬ 
sulting  technique,  referred  to  as  the  Gauss  Constrained  Technique,  was 
highly  efficient  and  required  very  few  objective  function  evaluations  to 
obtain  an  optimal  design.  Their  method  has  been  extended,  in  this  paper 
to  handle  nonlinear  constraints. 

Development  of  the  Method 

The  optimization  problem  consists  of  minimizing: 

f<;>  -  jT;  <i) 

where  $  is  a  vector  of  general  functions  of  the  variables  x.  Many  opti¬ 
mal  design  problems  have  objective  functions  of  this  form.  For  the 

derivation  of  the  optimization  technique,  the  0's  are  approximated  by 

linear  functions  of  x  (which  makes  F  quadratic)  of  the  form: 


where  J  and  c  are  a  constant  matrix  and  vector  respectively.  However 


the  technique,  once  derived,  will  be  applied  to  more  general  cases  where 
■+■  -*■ 

p  are  highly  nonlinear  functions  of  x.  The  constraints  are  approximated 
by  a  general  quadratic  of  the  form: 


gt(x) 


1 

2 


vr 

x  A 


-+T-+ 


.  x  +  b  .  x 
1  i 


d.  <  0 

l  — 


i=l . k 


(3) 


The  gradient  of  the  objective  function  (1)  is  given  bv: 


VF(x)  =  2J<£ 


(4) 


which  is  exact  whether  <f>  is  linear  or  not,  as  long  as  J  is  the  matrix 
of  first  partial  derivatives  of  <p .  The  matrix  of  second  partial  de¬ 
rivatives  of  the  objective  function  is: 


T 

G  =  2JJ 


(5) 


which  is  exact  only  when  <p  is  linear.  At  any  given  iteration  we  assume 
there  are  l  active  constraints  (ordered  such  that  j  =  1,2,..., £)  where 
l  <  k.  Thus,  the  Kuhn-Tucker  conditions  are: 

VF(x)  +  (A^x  +  b^)X^  +  ...  +  (A^x  +  =  0  (6) 

g j ( x )  =  -  x  A^x  +  bS  -  dj  =  0  j*l ,...,£  (7) 

Aj  >  0  j=l,...,£ 


Similar  to  the  development  of  Paradis  and  Willmert  [3],  the  new 


4 


method  is  generated  such  that  a  single  iteration  yields  the  optimal  de¬ 
sign  for  a  quadratic  objective  function  assuming  that  the  constraints 
active  at  the  optimal  design  are  also  active  at  the  starting  point.  If 

F(x)  is  quadratic  then  the  following  equation  is  valid  for  any  two  points 

•+•  , 
x  . ,  and  x  : 
u+1  v 

VF(x  ..)  -  VF(x  )  =  G(x  -  x  )  (9) 

V+l  V  V+l  V 


Solving  equation  (6)  for  VF(xv+^),  assuming  x^+^  is  the  optimal  design, 
and  substituting  into  equation  (9)  results  in  the  iterative  expression 
(after  solving  for  x^+^): 


Vi  ■ <G  +  -  ,r(\>1 


If  the  objective  function  has  the  special  form  of  equation  (1),  then 

the  G  matrix  in  the  iterative  equation  would  be  replaced  by  equation  (5), 
T 

i.e.  G  ■  2JJ  .  Otherwise,  the  method  could  still  be  applied  as  long  as 
the  matrix  of  second  partial  derivatives  G,  or  an  approximation  to  it, 
were  known. 

In  order  to  use  the  iterative  equation  (10)  the  unknown  vector  of 
Lagrange  multipliers  A  must  be  determined.  This  is  accomplished  by  sub¬ 
stituting  the  expression  for  xv+^>  equation  (10),  into  each  of  the  active 
constraints  of  equation  (7)  resulting  in  a  system  of  X,  nonlinear  equa- 
tions  in  i  unknowns  A  of  the  form: 


hM  ■{{((!<■  -  7F(Jv)])A1 

((G  +  jllAjXj)"llG*v  *  JjSjlj  '  VF(Jv)I) 

+  bj(G  +  )|1AjAj)'1[gJv  -  jIjSjXj  -  7F(JV)] 

-  dt  =  0  i«l,...,£  (11) 

This  system  of  equations  in  X  could  be  solved  using  several  dif¬ 
ferent  methods.  The  approach  used  in  this  research  was  as  follows.  At 
each  iteration  initially  the  Newton-Raphson  technique  was  applied.  If 
the  method  did  not  converge,  which  may  be  caused  by  the  fact  that  the 
equations  had  no  solution  or  that  the  method  simply  was  not  able  to 
locate  one,  then  another  approach  was  used.  In  this  case,  an  objective 
function  H  was  formed  which  was  the  sum  of  the  squares  of  the 
active  constraints  g^.  This  unconstrained  function  was  minimized  with 

respect  to  X.  In  this  work,  Powell's  method  was  used,  but  any  available 

*► 

technique  could  be  applied.  The  optimal  values  of  X  (whether  H  is  zero 
or  not  for  these  values  of  X)  were  then  used  in  the  iterative  equation 
(10). 

As  in  the  original  Gauss  Constrained  Technique  (as  well  as  the  gra¬ 
dient  projection  method)  the  criterion  for  dropping  a  constraint  from  the 
set  of  active  constraints  is  the  sign  of  x^-  At  any  iteration  the  con¬ 
straint  corresponding  to  the  most  negative  X  is  dropped  from  the  set  of 
active  constraints  and  a  new  vector  of  Lagrange  Multipliers  is  deter¬ 
mined.  This  procedure  is  repeated  until  all  Xi  >  0.  If  this  results  in 


no  active  constraints,  then  the  optimization  technique  continues  by  set- 
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ting  all  *  0  in  equation  (10).  In  this  case,  the  iterative  expression, 
can  be  shown  to  reduce  to  Gauss'  method  for  unconstrained  optimization: 


x  -  (JJT)-1J$ 
v 


(12) 


At  each  iteration  a  decision  must  be  made  as  to  whether  or  not  a  new 

constraint  should  be  added  to  the  set  of  active  constraints.  If  the 

direction  of  minimization  s  . ,  is  defined  as  the  direction  from  x  to 

v+1  v 

then  a  constraint  is  added  only  if  a  step  in  the  direction 

will  not  satisfy  the  constraint. 

Two  different  methods  were  used  to  handle  violated  constraints. 

The  first  approach  was  to  do  nothing  special  if  a  constraint  became 
violated  as  a  result  of  an  iteration  based  on  equation  (10).  Thus 
whether  a  constraint  is  only  active,  i.e.  *  0,  or  if  it  is  violated, 
i.e.  g.  >  0,  it  is  treated  the  same  in  the  optimization  technique.  The 
basis  behind  this  approach  was  the  assumption  that  eventually  the  method 
will  satisfy  all  constraints,  since  these  are  equations  (11)  which  the 
method  attempts  to  satisfy  at  each  iteration.  This  approach  allows  the 
optimization  to  start  at  any  design  whether  it  satisfies  the  constraints 
or  not. 

The  other  method  of  handling  violated  constraints  was  to  start  at 

a  point  which  satisfied  all  constraints.  Then  if  a  step  is  taken  using 

equation  (10)  such  that  gr(x^+^)  >  0,  a  new  design  x^+^  was  found  along 

the  line  connecting  x  and  x  , ,  such  that  g  (x  .,)  =  0.  This  was 

v  v+1  °r  v+1 

accomplished  by  stepping  back  from  a  violated  constraint  using  the 
approximate  equation: 


L> 


(x  -  X  )g  (x  ) 
-*■  v+1  v  6r  V 

X  •  '■ 

g  (x  1 )  -  g  (x  ) 
°r  v+1  ®r  v 


The  design  x^+^  was  then  used  in  the  next  iteration. 

Although  the  technique  was  developed  assuming  linear  functions  <J>(x), 
i.e.  F(x)  quadratic,  the  method  is  applicable  to  problems  where  the  $s 
are  general  functions  of  x.  An  especially  important  characteristic  for 
mechanism  design  is  that  the  technique  requires  a  total  of  only  one  more 
objective  function  evaluation  than  iterations  to  obtain  an  optimal  design, 
unlike  the  GRG  Technique  of  Lasdon  et.al  [2],  where  a  step  length  deter¬ 
mination  is  required  at  each  iteration. 


Example  Problems 

Several  examples  were  considered  in  this  work,  two  of  which  are 
presented  here. 

The  first  problem,  which  is  a  modified  Rosen-Suzuki  test  problem, 
consists  of  a  quadratic  objective  function  in  eight  variables  with 
seven  quadratic  inequality  constraints  as  follows: 


Minimize 


2  3212,-2  2  52  2.32.. 

F(x)  *  x^  +  — x^  +  -^x^  +  3x^  +  2x^  +  -^Xg  +  x^  +  -^Xg  +  2x^^ 


3x2  +  5x4  +  x5  +  ?x6  -  2xy  -  xg 


4 


Tables  1  and  2  show  the  results  of  the  optimization  from  two  dif- 
ferent  starting  points.  Three  different  techniques  were  applied.  The 
first  two  were  variations  of  the  Gauss  Nonlinear ly  Constrained  method. 

The  first  one,  labeled  GNLC.NS  in  the  tables,  ignored  violated  con¬ 
straints,  while  the  second,  labeled  GNLC,  treated  violated  constraints 
using  equation  (13).  Both  of  these  methods  were  compared  with  the  Gen¬ 
eralized  Reduced  Gradient  technique,  labeled  GRG. 

Table  1  shows  the  results  from  a  starting  point  where  no  constraints 
were  active  or  violated.  All  three  methods  produced  the  same  optimal 
design.  The  number  of  iterations  NI,  and  number  of  function  evaluations 
NF  for  the  GNLC  methods  were  considerably  less  than  for  the  GRG.  Table  2 
shows  the  results  from  a  starting  point  in  which  three  constraints  were 
violated,  while  one  was  active.  For  this  starting  point  only  the 
GNLC.NS  technique  was  applicable.  Again  the  number  of  iterations  and 
function  evaluations  were  less  than  for  the  GRG. 

The  second  example  consists  of  the  design  of  a  four-bar  mechanism 
such  that  the  coupler  point  (XX, YY)  of  the  mechanism  closely  generated 
the  curve  defined  by  the  eight  points  as  shown  in  Figure  1.  The  ob¬ 
jective  function  to  be  minimized  is  the  sura  of  the  distances  (squared) 
between  the  desired  curve  and  the  actual  curve  generated  by  the  mechanism 

-*■  ®  2  2 
F(x)  =  .^.(XX,  -  XG.)  +  (YY.  -  YG  ) 
i*l  i  l  i  i 


XXi  *  x2C08Yi  +  XgCOS!^  +  x10 


where  s 


tan 


-1 


x-sin(Y<  -  x.) 


Xj^  -  X2C08(Yi  -  XgJ 


cos 


-1 


_2  ^  2  ^  2  2-  ,  % 
Xj^  +  X2  +  x^  “  x^  -  2x^X2COS(y^  "  Xg) 


r  2  2 

2x3  xi  +  x2  ‘  2x1x2coa(Yi  “  xfe) 


1/2 


Yi  -  xg  +  Ay£ 

and  (XG,YG)  are  the  coordinates  of  the  desired  curve. 

The  design  variables  x^,...,x^0  are  link  lengths,  orientation  of 
some  of  the  links  as  well  as  the  location  of  the  crank  pin.  Note,  the 
objective  function  is  a  highly  nonlinear  function  of  the  design  vari¬ 
ables.  This  is  the  same  problem  that  Paradis  and  Willmert  [3]  solved, 
however,  two  nonlinear  constraints,  g^3  and  g^,  have  been  introduced  here 
to  constrain  the  location  of  the  crank  pin  to  two  user  defined  circular 
regions.  Thus,  our  design  problem  contains  ten  variables,  x^,...,x^q, 
and  fourteen  inequality  constraints  g^(x)  <  0  i*l,...,14.  The  con¬ 


straints  are: 


g^x)  ■  1  -  x2  <  0 


g2<x) 

83(5t) 

84<*> 

g5(£) 

ha) 

g7(^) 


1  -  Xj  <  0 

x2  -  xl  <  0 

x2  “  x3  —  0 
x2  x4  —  0 
x3  +  x2  -  x3 
x2  +  x3  -  xx 


x.  <  0 
4  — 


x,  <  0 
4  — 


gg(5r)  »  x2  +  x4  -  x1  -  x3  <  0 
g^(^)  *  x^  -  30  <  0 


g1Q(£)  *  x3  -  30  <  0 

gn(*>  -  x4  -  30  <  0 

g12^)  *  x$  -  30  <  0 

g13(*)  »  (xiq  -  CXi>2  +  (x9  -  CY1)2 

-  R32  <  0 

g14(x)  »  (x1Q  -  CX2)2  +  (x9  -  CY2>2 

-  r22  <  0 


where  CX^,  CY^,  ,  CX^,  CY2,  and  R2  were  4,  4,  3,  7,  7,  and  3  re¬ 
spectively. 

The  resulting  optimal  designs  using  the  GNLC  and  GRG  methods  are 
shown  in  Table  3.  The  optimal  mechanism  is  shown  in  Figure  2.  Again 
the  Gauss  Nonlinearly  Constrained  method  required  significantly  fewer 
objective  function  evaluations. 


Conclusions 

The  Gauss  Nonlinearly  Constrained  method  is  an  effective  technique 


of  solving  nonlinear  design  problems.  It  is  particularly  efficient 
for  cases  in  which  the  evaluation  of  the  objective  function  is  very 
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time  consuming.  Although  in  some  instances  considerable  calculations 
must  be  done  per  iteration,  the  amount  required  is  still  insignificant 
compared  to  that  required  to  do  just  one  analysis  of  a  complex  mechanical 
or  structural  system. 

The  method  has  been  applied  to  a  variety  of  problems  consisting  of 
objective  functions  of  various  complexities.  In  all  cases  it  has  worked 
well.  However,  it  has  received  only  limited  application  to  problems 
involving  constraints  that  are  more  than  quadratically  nonlinear.  It 
appears  that  further  research  is  required  for  such  cases. 
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Table  1:  Results  of  Rosen-Suzuki  Test  Problem 
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Table  2:  Results  of  Rosen-Suzuki  Test  Problem 
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ABSTRACT 

A  finite  element  approach  is  presented  in  this  paper  for  the  nonlinear 
vibrational  analysis  of  planar  mechanisms.  The  analysis  takes  into  account 
the  effects  of  material  and  geometric  nonlinearities  on  the  dynamic  behavior. 
The  geometric  nonlinearities  included  in  this  study  are  due  to  stretching  of 
the  neutral  axis  and  the  curvature-displacement  nonlinearity,  both  caused  by 
large  deformations.  The  material  nonlinearity  is  due  to  a  nonlinear  stress- 
strain  relationship  of  the  Ramberg-Osgood  type.  The  analysis  presented  here 
makes  use  of  hermite  polynomials  which  ensure  compatibility  of  curvature 
between  elements.  Using  a  variable  correlation  table,  a  global  system  of 
nonlinear  equations  are  derived  in  terms  of  the  global  unknowns  and  the 
kinematics  of  the  mechanism.  A  harmonic  series  technique  is  then  used  to 
obtain  the  steady  state  solutions  to  this  system  of  nonlinear  equations. 
Numerical  results  are  presented  for  an  example  mechanism  and  the  effects  of 
the  nonlinearities  are  discussed. 


INTRODUCTION 


The  importance  of  flexibility  of  linkages  on  the  performance  of  high¬ 
speed  minimum-mass  mechanisms  is  well  recognized.  A  considerable  amount  of 
research  has  been  done  in  this  area  in  the  last  two  decades.  While  it  is 
desirable  to  develop  analytical  and  numerial  procedures  that  enable  the 
design  of  rigid  link  mechanisms  and  robots  to  perform  a  given  function  with 
specified  reliability,  it  is  also  important  to  evaluate  the  effects  of 
flexibility  of  elastic  members  on  their  performance.  It  is  known  that  a 
mechanism  designed  for  operation  at  low  speeds  may  not  perform  satisfactorily 
at  high  speeds  due  to  the  effects  of  large  inertia  forces  and  resulting 
elastic  deformations.  Thus  it  becomes  necessary  to  include  in  the  dynamic 
analysis  of  mechanisms,  not  only  the  effect  of  the  rigid  body  motion,  but 
also  the  flexibility  of  the  linkages. 

Most  of  the  previous  investigations  in  the  area  of  elastic  analysis  of 
mechanisms  have  been  carried  out  within  the  framework  of  the  linear  theory 
[1-16].  However,  Viscomi  and  Ayre  [17]  used  a  Galerkin-type  nonlinear 
analysis  procedure  to  study  the  vibrations  of  a  slider-crank  mechanism.  A 
later  work  by  Sadler  and  Sandor  [18]  used  the  lumped  parameter  approach  to 
a  nonlinear  dynamic  model  of  an  elastic  linkage.  The  mechanism  analysed  in 
this  paper  was  a  general  four-bar  linkage  and  the  analytical  model  inc’uded 
the  response  coupling  associated  with  both  tne  transmission  of  forces  at  the 
pin  joints  and  the  dependence  of  the  undeformed  motion  of  a  link  on  the 
elastic  motion  of  other  links.  A  finite  element  analysis,  w<  th  the  aid  of 


I  *>  *  %  *- >  L". 
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Che  piecewise  linear  method  of  Martin,  was  used  by  Sevak  and  McLarnan  [19] 
to  carry  out  the  nonlinear  analysis  of  a  mechanism.  Further  nonlinear  work 
dealing  with  the  vibrations  of  elastic  mechanisms  are  reported  in  References 
[20-22],  In  a  recent  investigation,  Thompson  and  Sung  [23]  used  a  varia¬ 
tional  formulation  for  the  nonlinear  finite  element  analysis  of  planar 
mechanisms  considering  geometric  nonlinearities.  Some  experimental  results 
were  presented. 

This  paper  is  concerned  with  the  nonlinear  vibrational  analysis  of 
general  planar  mechanisms.  A  finite  element  method  is  used  which  includes 
the  effects  of  both  geometric  and  material  nonlinearites .  The  geometric 
nonlinearities  included  in  this  study  are  due  to  stretching  of  the  neutral 
axis  with  partially  constrained  ends  and  a  general  curvature-displacement 
relationship,  both  caused  by  large  deformations.  The  material  nonlinear¬ 
ity  is  of  the  Ramberg-Osgood  type  with  three  parameters  to  represent  the 
nonlinear  stress-strain  relationship  [24-26].  Additional  effects  consid¬ 
ered  are  transverse  shear  and  rotatory  inertia  and  changes  in  cross-section 
due  to  realistically  proportioned  members.  The  governing  nonlinear  differ¬ 
ential'  equations  are  derived  for  each  element  in  terms  of  the  axial  and 
transverse  deformations,  rotations,  curvatures,  and  shear  deformation  angles. 
These  equations  are  then  assembled  with  the  aid  of  a  variable  correlation 
table  and  the  resulting  global  system  of  equations  is  solved  using  an 
iterative  technique  based  on  a  harmonic  series  solution  procedure. 
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FINITE  ELEMENT  FORMULATION 


A  finite  element  method  is  presented  below  for  the  nonlinear  analysis 
of  a  general  closed  looped  mechanism.  The  mechanism  can  be  composed  of 
various  combinations  of  simple  four  bar  chains,  frame  elements,  sliders 
moving  on  fixed  references,  or  sliders  moving  on  rotating  links.  Each  link 
is  divided  into  one  or  more  elements  with  each  element  having  the  local 
coordinate  system  as  shown  in  Figure  1.  If  a  slider  is  present,  the  masses 
M^  and  M^  are  located  at  ends  1  and  2,  as  shown.  The  length  of  the  element 
A  is  constant  except  for  links  with  sliders  moving  along  them. 

The  displacement  vector  of  any  point  (a)  on  the  element's  neutral  axis 
is  given  bys 

S  ■  (X^  cosy  +  Y^  siny  +  x  +  u)i  +  (Y^  cosy  -  siny  +  w)j  (1) 

where  X^  and  Y^  are  the  coordinates  of  end  1  of  the  element  given  by  the 
rigid  body  motion.  The  coordinate  x  is  measured  along  the  element's  neutral 
axis  from  1  to  2  and  y  is  the  angle  between  the  rigid  body  position  and  the 
X-axis.  The  axial  and  transverse  displacements  of  point  (a)  from  the  rigid 
body  position  are  given  by  u  and  w,  respectively.  This  equation  takes  into 
account  both  the  rigid  body  motion  and  the  elastic  displacements  and  defines 
the  position  of  any  point  along  the  neutral  axis. 

Differentiating  Equation  1  with  respect  to  time  yields  the  velocity  of 
any  point  (a).  The  unit  vectors  i  and  j  move  with  the  coordinate  system  and 


vary  with  time.  The  angular  velocity  of  any  differential  line  segment  on 
the  neutral  axis  of  the  element  is  given  by: 


+  w 


*xt 


(2) 


where -Y>t  is  the  derivative  of  y  with  respect  to  time  and  w>xt  Is  the  de¬ 
rivative  of  the  transverse  displacement  with  respect  to  the  local  coordinate 
x  and  time  t. 

The  kinetic  energy  due  to  rotation  of  the  element  is  given  by: 


K.E 


‘R 


^  1  2  2 

J  p Ax I s > t !  +  pI2(y,t  +  w,xt}  dx  dy 


+lM1|S,tl2 


+  k  Mjs,  I 

x=o  2  2 1  ’  t 1 


2 

x“A 


(3) 


where  p  is  the  mass  density  and  and  Iz  are  the  cross  sectional  area  and 

,  2 
moment  of  inertia  of  the  element  respectively.  The  term  plz(y,t  +  w»xt^  /2 

represents  the  effect  of  rotatory  inertia.  The  kinetic  energy  due  to  beam 

bending  associated  with  transverse  shear  is  [27]. 

h 

1  ^  A  2 

K.E.b  =  2  J  J  {a,tIz>dx  dy  (4) 

ho 

"  2 

where  a  is  a  measure  of  the  transverse  shear  angle. 
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u  -  J  fa  de  dvol  +  J  |  t 


vol  o 


xy  xy 


where  o,  e,  t  ,  and  y  are  the  normal  stress,  normal  strain,  shear  stress 
xy’  'xy 

and  shear  strain,  respectively.  For  a  nonlinear  material  of  the  Ramberg- 
Osgood  type  [24-26],  the  relationship  between  stress  and  strain  is: 


o  «  A  e  -Be 


where  A  corresponds  to  Young's  modulus  E,  and  Bem  represents  the  nonlinear 
term.  A,  B  and  m  are  constants  for  the  particular  material  being  considered 
The  above  relationship.  Equation  6,  is  valid  only  for  positive  strains.  If 
the  strain  is  negative,  the  following  expression  is  used: 


aaAe+B(-e)  if  e  <  0 


The  change  in  sign  of  the  nonlinear  term  results  in  the  same  overall  effect 
on  the  stress-strain  relationship  as  for  positive  strain,  i.e.  either  hard¬ 
ening  or  softening  depending  on  the  values  of  B  and  m. 

Using  the  Ramberg-Osgood  relationship,  the  following  expression  for  the 
strain  energy  is  obtained  for  positive  strain  e  >  0: 


h  h 

2  2  ^ 

u  =  J  /A4  A  £2  -  B  e1"*1  dx  dy  +  T  \  \  G  y2  dx  dy  (8) 


where  x  is  the  axial  coordinate,  y  is  the  transverse  coordinate,  and  G 

3  xy 

is  the  shear  modulus.  When  e  <  0  the  equation  is: 


2  A 

U  *  J  J  {-^  A  e2  -  B  (-d^ldx  dy 
h  o 
"  2 

h 

2  A1  2 

+  J  J  4  G  y  dx  dy  (9) 

1  ;  2  xy  > 
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The  nonlinear  expression  for  the  curvature  R  of  a  planar  static  beam 
undergoing  large  deformations  is: 


1 

R 


w. 


xx 


M  ,  2,3/2 

(1  +  w,  ) 

X 


(10) 


Thus  the  strain  is: 


y  w, 


xx 


(1  +  w 


2,3/2 

’x' 


(11) 


Combining  the  geometric  nonlinearities  due  to  stretching  of  the  neutral  axis 
and  the  curvature-dispalcement  nonlinearity,  results  in  the  expressions  for 
normal  and  shear  strain: 


u,  +  4  D 
x  2  s 


y  a. 


w, 


( 1  +  D,  w 


2}3/2 


(12) 


Substituting  these  expressions  for  strain  into  the  strain  energy 
tions  (8)  and  (9)  produces,  for  e  >_  0 


S.E. 


J  J  b{|  A  (u,v  +  \  w,v  - 


y  cti 


h  o 
2 


x  ^  „  2.3/2 

(1  +  Dl  w,  ) 

D  A 


Dk  1  2  y  “’x 

+  B  (u,  +  ±  D  v,  - 

m+1  '  x  2  s  x 


(1  +  D.  w,2)3/2 
b  x 


m+1 


+  I  G  (w,  +  a)2}dx  dy 

2  xy  x  J 


and  when  z  <  0 


S.E. 


-  J  I  b  {?  A  (u,y  +  |  Dg  w,x 
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h  o 
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S  X 
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+  4  G  (w,  +  ct)^}dx  dy 

2  xy  x  J 


Equa 


(14) 


(15) 


In  Equations  (12),  (14)  and  (15),  D  ,  D,  and  D  are  tracing  constants 


representing  the  effects  of  material  nonlinearity,  geometric  nonlinearity 
due  to  curvature,  and  geometric  nonlinearity  due  to  stretching  of  neutral 
axis,  respectively.  Each  tracing  constant  is  equated  to  unity  when  that 
particular  nonlinearity  is  being  considered  and  is  equated  to  zero  when  it 
is  not. 

In  order  to  represent  realistically  proportioned  members,  changes  in 
cross  section  are  included.  Each  element  is  divided  into  sections  of  varying 

m 

lengths  with  constant  area.  The  intergrations  involved  in  the  element  equa¬ 
tions  are  carried  out  in  a  piecewise  fashion  with  the  area  in  each  section 
taken  as  a  constant.  This  procedure  provides  a  reasonable  approximation  of 
variable  cross  sectional  members  without  having  to  resort  to  large  numbers 
of  elements. 

The  Lagrange  function  L  is  defined  as: 

s  \ 

L  ■  kil  1-E1  <k-e-r  +  k-e-b  ■  s-E->ik  (16) 


where  N,  is  the  total  number  of  elements  in  the  kth  link  and  S  is  the  num- 
k 

ber  of  links  in  the  mechanism.  Substituting  Equations  (3),  (4)  and  (14) 
into  Equation  (16),  the  Lagrangian  L  can  be  expressed  in  terms  of  the  dis¬ 
placements  u,  w,  the  shear  angle  a,  and  the  rigid  body  motion. 

Hermite  polynomials  are  used  to  approximate  u,  w,  and  a  in  order  to 
satisfy  the  boundary  conditions  of  various  types  of  mechanisms  easily  and 
to  ensure  interelement  compatibility.  The  axial  deformation  u  is  approxi¬ 
mated  by  a  linear  shape  function  given  by 

I 

u  -  Vc)  Nl^>  +  VC)  N2(x)  (17) 

I 
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Similarly,  fifth  degree  polynomial  shape  functions  are  used  to  approximate 
the  transverse  deformation  w: 


w  =  WjCt)  +  6^(0  H21(x)  +  m^t)  H31(x) 

+  W2(t)  H12(x)  +  §  (t)  H22(x)  +  m2(t)  H32(x)  (18) 

where  6  and  m  are  continuous  between  elements.  The  shear  angle  a  is  also 
approximated  by  a  fifth  degree  polynomial  in  order  to  make  it  compatible 
with  the  transverse  displacement  w.  Therefore  a  is  assumed  to  be: 

a  *  a^(t)  H^x)  +  iL(t)  H^Cx)  +  X^(t)  H^x) 

+  i2<t)  H12(x)  +  J2(t)  H22(x)  +  X  (t)  H32(x)  (19) 

where  \p  and  X  are  the  first  and  second  derivatives  of  a,  respectively. 

The  Hermite  polynomials  are  given  by: 

N^(x)  *  1  -  e 

N2(x)  *  e  (20) 

Hu(x)  -  1  -  10e3  +  15e4  -  6e5 
H2i(x)  *  A(e  -  6e3  +  8e4  -  3e3) 

H3i(x)  *  -  3e3  +  3e4  -  e3)/2 

H  (x)  *  10e3  -  15e4  +  6e5 

H22^X^  *  A(_4e^  +  7e4  -  3e3) 

«32(x)  -  A2(e3  -  2e4  +  e5)/2  (21) 
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where,  e  =  x/A. 


A  transformation  of  coordinates  is  now  introduced  to  change  from  the 
moving  coordinate  system  associated  with  the  elements  to  global  coordinates. 
Only  U^,  Wj,  U2  and  W2  need  to  be  transformed.  The  other  coordinates  are 
angles  or  derivatives  of  angles  which  are  not  directional  on  the  X,  Y  co¬ 
ordinate  system  used.  The  transformations  are: 


cos<^  -  sin<^ 

W1  *  U1  sin$i  +  cos<)>^ 

U,  =  U2  cos<j>2  “  W2  sin^2 

W2  =  U2  sin«j>2  +  co8$2  (22) 


For  pin  connections  the  transformation  angles  and  are  set  equal  to  -y 
(the  rigid  body  angle)  which  transforms  the  coordinates  back  to  the  global 
coordinates.  For  sliders  moving  on  rotating  links  the  transformation  becomes 
more  involved.  In  this  case  the  deformation  of  the  driver  link  must  be 
transformed  to  correspond  to  the  axial  and  transverse  deformation  of  the 
rotating  link  [14]. 

Substituting  the  expressions  form  Equation  (22)  into  Equations  (17)  and 
(18),  the  global  coordinates  for  the  system  are  then: 


q  =  [Uj  ex  mx  a1  A ^  U2  W£  $2  m2  a2  * 2  X2^ 


(23) 


The  Lagrangian  function  is  then  written  in  terms  of  the  transformed  element 
coordinates.  Differentiating  the  Lagrangian  with  respect  to  the  element  co¬ 
ordinates,  the  following  element  equations  are  obtained: 


3 


E 


2' 


dt  3q,t  3q 

In  differentiating  the  expressions  for  the  kinetic  and  strain  energies  in 
Equations  (3),  (4)  and  (14)  it  must  be  kept  in  mind  that  A  which  is  the 
upper  limit  of  integration  is  a  function  of  time.  The  operations  carried 
out  in  Equation  (24)  results  in  a  system  of  nonlinear  element  differential 
equations.  Assembling  the  element  matrices  for  the  particular  mechanism 
being  solved  results  in  the  global  system  of  equations: 

M  Q,  +CQ,.  +  (K  +  K  )  Q  =  F(t)  (25) 

c  t.  c  n 

The  M,  C,  K  and  K  matrices  are  all  functions  of  time.  The  C  matrix  re- 
e  n 

suits  from  the  kinetic  energy  of  the  system.  No  damping  was  included  in 

the  formulation  of  the  problem.  The  C  Q1(.  term  was  found  to  be  small  and 

thus  was  ignored  in  the  analysis.  The  matrix  Ke  is  the  linear  portion  of 

the  total  stiffness  matrix.  It  is  a  function  of  rigid  body  motion  but  not 

a  function  of  the  deformations  Q.  The  matrix  Kn>  however,  is  the  nonlinear 

portion  of  the  stiffness  matrix.  It  is  a  function  of  the  deformations. 

Equation  (25)  is  thus  a  nonlinear  system  of  equations. 

The  derivation  of  the  finite  element  Equation  (25)  is  based  on  the 

assumption  of  positive  strains  e.  If  the  strain  is  negative  a  similar 

derivation  is  possible,  based  on  Equation  (15)  for  the  strain  energy  rather 

than  Equation  (14).  The  only  difference  in  the  resulting  Equation  (25)  is 

in  the  stiffness  matrix.  Wherever  an  em  ^  occurs  in  the  original  formula- 

m-1 


tion,  it  becomes  (-e) 


for  negative  strains.  All  other  negative  signs 


resulting  from  the  introduction  of  -e  cancelled  out  in  the  differentiations 
required.  Thus  in  order  to  handle  both  positive  and  negative  strains  the 
terms  involving  em  ^  in  the  stiffness  matrix  were  replaced  by  lef1. 

In  order  to  solve  the  nonlinear  system  of  Equation  (25)  an  iterative 
approach  was  used.  First  the  equations  were  solved  using  the  linear  terms 
only,  i.e.  the  matrix  was  ignored.  This  was  accomplished  by  setting 
all  of  the  tracing  constants  D^,  and  equal  to  zero.  The  solution  Q 
for  the  linear  equations  was  then  used  to  determine  values  for  the  nonlinear 
stiffness  matrix  K^.  Equation  (25)  was  then  solved  again  for  new  values 
for  Q,  and  the  process  repeated.  Experience  showed  that  this  procedure 
converged  in  from  3  to  5  iterations.  To  solve  Equation  (25)  for  Q(t),  for 
particular  K^,  a  harmonic  series  solution  method  was  used  similar  to  that 
of  Bahgat  and  Willmert  (14J.  This  approach  overcomes  problems  with  sta¬ 
bility,  due  to  the  time  varying  nature  of  the  matrices,  that  sometimes  re¬ 
sult  from  an  eigenvalue  technique.  The  steady  state  solution  is  obtained 
without  adding  artificial  damping.  The  solution,  without  the  C  Q,t  term, 
is  given  by: 


N  22-1- 

Q(t)  =  Z  (K  -  n  w  M)  (A  cos  nut  +  B 
n=o  n  r 


sin  n^l 


where w  is  the  input  crank  speed,  and  A  and  B  are  solutions  to  the  linear 

n  n 

equations : 


N  N-  1  _ 

F(t,  )  =  Z  A  cos  nujt,  +  Z  B  sin  nut, 
k  nHo  n  k  n=o  n  1 


for  k  =  0,1,  .. . , 2N- 1 
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where  Aq  is  set  equal  to  zero.  The  values  of  t^  are  the  times  at  2N  equal 
time  increments  per  revolution  of  the  input  crank  given  by 

ITlc 

tfc  -  —  for  k  =  0,1,..., 2N-1  (28) 

Computational  experience  indicates  that  a  fairly  accurate  solution  is  ob¬ 
tained  using  only  a  few  terms  in  Equation  (26).  As  the  number  of  terms  in- 

2  2 

creases  the  components  of  the  matrix  (K  -  n  w  M)  grow  and  thus  the  in- 
2  2  - 1 

verse  (K  -  n  u  M)  becomes  small.  The  summation  can  therefore  be 
truncated  to  reduce  computational  time. 

The  stress  in  the  links  is  calculated  by  evaluating  the  strains  from 
Equation  (12).  The  stress  can  then  be  determined  at  any  point  in  an  element 
using  Equations  (6)  or  (7).  To  find  the  maximum  stress  in  an  element  the 
maximum  strain  must  be  found.  Setting  the  first  derivative  of  Equation  (12) 
to  zero  and  solving  the  resulting  expression,  the  position  of  the  maximum 
strain  is  determined.  Once  the  location  is  known,  the  maximum  strain  and 
stress  can  be  evaluated. 

The  above  formulation  is  based  on  the  use  of  the  shear  angle  a,  which 
is  appropriate  particularly  for  short  members.  For  long  slender  links  this 
quantity  is  not  required.  The  elimination  of  a  reduces  the  size  of  the 
problem  considerably  since  the  nodal  deformations  a^,  iji  ,  X^,  and 

A2  would  no  longer  be  present.  For  long  slender  members: 


-  w,. 


(29) 
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.  **  ■ " .  •  fc-*  ■-*  o 


Using  this  expression,  the  equation  for  strain  energy  (14)  for  positive 
shear  e  reduces  to: 


S.E.  * 


h 

2 

J 

_h 

"2 


A  (u. 


2  °s 


2 

’x 


y  w.. 


(1  +  D,  w 


,V/2 


nH-1 


(u,  +  \  D  w,2  + 

x  2  s  x 


y  w». 


(1  +  D  w,2)3/2 
b  x 


)m+1}dx  dy 


(30) 


A  similar  expression  exists  for  negative  strain.  The  kinetic  energy  also 
changes  if  a  is  not  present.  The  energy  associated  with  transverse  shear, 
Equation  (4),  is  eliminated  and  thus  Equation  (3)  represents  the  total 
kinetic  energy  of  the  element.  Using  a  procedure  similar  to  the  method 
outlined  above,  vibrational  equations  of  the  same  form  as  Equation  (25)  can 
be  obtained,  but  they  will  be  smaller  in  size.  However,  nonlinear  terms 
still  exist  in  the  stiffness  matrix  due  to  the  material  and  geometric  non- 
linearities  The  method  of  solution  is  thus  identical  to  that  outlined 
earlier . 
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EXAMPLE  PROBLEM 


The  following  example  is  presented  to  illustrate  the  method  of  solution. 
The  nonlinearities  due  to  neutral  axis  stretching,  curvature-displacement 
and  stress-strain  relationship  are  all  considered.  A  four-bar  linkage, 
as  shown  in  Figure  2,  is  used  as  the  example  with  all  the  members  flexible 
and  made  of  the  same  material.  The  data  for  the  mechanism  is 

Length  of  input  crank  (AB)  »  5.0  in 
Length  of  coupler  (BC)  *  11.0  in 
Length  of  rocker  (CD)  *  10.5  in 
Fixed  distance  (AD)  ■  10.0  in 
Cross  section  of  links  ■  rectangular 
Height  of  rectangle  *  1.0  in 
Width  of  rectangle  *  0.25  in 

The  initial  position  of  the  input  crank  is  zero  degrees  at  t  *  0  and  the 
direction  of  rotation  is  counterclockwise.  The  mechanism  is  divided  into 
three  elements  with  each  link  in  the  mechanism  taken  as  an  element.  The 
boundary  conditions  are  that  only  moment  and  shear  terms  exist  for  the  input 
crank's  driven  end  (A).  For  the  pin  connections  between  links,  there  are 
deformations,  rotations  and  shear  terms,  and  for  the  rocker's  fixed  point 
there  are  only  rotation  and  shear  terms. 

First,  the  deformations  in  the  mechanism  were  determined  with  the  shear 
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angle  a  present.  In  this  case  the  crank  link  was  rotated  at  100  rad/sec. 
The  material  properties,  approximating  aluminum,  were  as  follows: 

A  -  10.87xl06  lb/in2 
B  »  0.8387X1011  lb/in2 
m  *  3.0 

Mass  density  =  0.0002536  lb-sec2/in^ 

Three  separate  procedures  were  used  to  obtain  numerical  results.  First 
the  problem  was  solved  using  the  linear  analysis  method  of  Bahgat  and 
Willmert  [14],  with  E  =  A.  Next  the  method  of  this  paper  was  used  with  the 
tracing  constants  equal  to  2ero.  Thus  a  linear  analysis  was  obtained. 
Finally  the  method  was  applied  with  all  tracing  constants  equal  to  one, 
i.e.  a  full  nonlinear  analysis.  A  representative  deformation  as  a 
function  of  crank  position  is  shown  in  Figure  3.  This  is  the  horizontal 
deformation  of  the  free  end  of  the  crank  link.  As  can  be  seen,  the  three 
curves  are  very  similar.  The  effect  of  the  shear  angle  a  is  to  increase 
the  deformation  slightly.  For  this  slow  speed  the  linear  and  nonlinear 
analyses  were  almost  identical. 

The  same  problem  was  solved  again  at  a  higher  speed  of  200  rad/sec. 

The  resulting  deformation  is  shown  in  Figure  4.  As  can  be  seen,  high 
frequency  oscillations  started  to  appear,  with  greater  separation  between 
the  three  analyses.  At  even  higher  speed  these  oscillations  became  more 
predominant  to  the  point  of  instabilities  in  the  analysis  at  very  high 
speeds . 
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The  revised  form  of  Che  analysis  equations  was  considered  next,  i.e. 

Che  form  without  the  strain  angle  a.  Here  a  crank  speed  of  150  rad/sec  was 
used.  A  comparison  was  made  of  the  effects  of  the  various  nonlinearities 
on  the  deformations  and  stresses  as  compared  to  the  linear  analysis. 

Figures  5  and  6  show  a  comparison  of  the  linear  and  nonlinear  deformations 
U2  (the  horizontal  deformation  of  the  free  end  of  the  output  link)  caused 
separately  by  geometric  and  material  nonlinearities.  Figures  7  and  8  show 
the  maximum  stresses  in  the  connecting  link  of  the  mechanism.  As  expected, 
the  material  nonlinearity  of  the  Ramberg-Osgood  type  results  in  deformations 
which  are  greater  in  magnitude  than  those  obtained  using  a  linear  elastic 
model.  The  maximum  stress  decreased  due  to  the  presence  of  the  term  Be 111 
subtracted  from  the  linear  stress  expression. 

The  geometric  nonlinearities  considered,  namely  curvature  displacement 
and  stretching  of  the  neutral  axis,  both  due  to  large  deformations  produced 
mixed  results  with  deformations  reduced  at  some  points  and  increased  at  other 
points.  The  effect  of  the  geometric  nonlinearities  would  be  expected  to 
produce  a  stiffening  of  the  members  [27]  of  the  mechanism  and  thus  produce 
smaller  deformations.  The  increased  deformations  in  this  case  might  be  due 
to  the  fact  that  the  deformations  are  in  relationship  to  the  entire  mechanism 
and  not  just  to  the  individual  beam  element. 


CONCLUSIONS 


The  nonlinear  analysis  procedure,  using  a  finite  element  technique,  is 
an  effective  method  of  calculating  the  steady  state  deformations  and  stresses 

in  a  mechanism.  Significant  differences  can  occur  between  the  linear  and 
nonlinear  approaches.  This  was  particularly  true  for  the  stresses  in  the 
example  considered  in  this  work.  Research  is  still  needed  on  the  overall 
effect  of  the  shear  angle  a,  and  a  more  complete  picture  of  the  nonlinear 
terms  in  the  analysis  would  be  of  value.  Additional  nonlinear  effect  should 
also  be  investigated,  such  as  separation  of  the  translations  of  one  link 
due  to  large  deformations  of  the  other  links. 
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Figure  3  -  Deformations  U.  for  Aluminum  at  100  rad/sec 


Figure  5  -  Comparison  of  Deformations  at  150  rad/sec  for  Geometric  Nonlinearities 


igure  6  -  Comparison  of  Deformations  U„  at  150  rad/sec  for  Material  Nonlinearit 
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Figure  8  -  Maximum  Stresses  in  Element  2  at  150  rad/sec  for  Material  Nonlinearity 


